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Abstract 



In [C. W. Gear, T. J. Kaper, I. G. Kevrekidis, and A. Zagaris, Projecting to 
^ i a Slow Manifold: Singularly Perturbed Systems and Legacy Codes, SIAM J. 

Appl. Dyn. Syst. 4 (2005) 711-732], we developed a class of iterative algorithms 
■ within the context of equation-free methods to approximate low-dimensional, 

attracting, slow manifolds in systems of differential equations with multiple 
\ time scales. For user-specified values of a finite number of the observables, the 

' m— th member of the class of algorithms (m = 0, 1, . . .) finds iteratively an 

. approximation of the appropriate zero of the (m + 1)— st time derivative of the 

remaining variables and uses this root to approximate the location of the point 
on the slow manifold corresponding to these values of the observables. This 
article is the first of two articles in which the accuracy and convergence of the 
5^ \ iterative algorithms are analyzed. Here, we work directly with explicit fast-slow 

systems, in which there is an explicit small parameter, e, measuring the separa- 
tion of time scales. We show that, for each m = 0, 1, . . ., the fixed point of the 
iterative algorithm approximates the slow manifold up to and including terms of 
C(e™'). Moreover, for each m, we identify explicitly the conditions under which 
the m— th iterative algorithm converges to this fixed point. Finally, we show 
that when the iteration is unstable (or converges slowly) it may be stabilized 
(or its convergence may be accelerated) by application of the Recursive Projec- 
tion Method. Alternatively, the Newton-Krylov Generalized Minimal Residual 
Method may be used. In the subsequent article, we will consider the accuracy 
and convergence of the iterative algorithms for a broader class of systems — in 
which there need not be an explicit small parameter — to which the algorithms 
also apply. 
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1 Introduction 



The long-term djTiamics of many complex chemical, physical, and biological systems 
simplify when a low- dimensional, attracting, invariant slow manifold is present. Such 
a slow manifold attracts all nearby initial data exponentially, and the reduced dy- 
namics on it govern the long term evolution of the full system. More specifically, a 
slow manifold is parametrized by observables which are typically slow variables or 
functions of variables. All nearby system trajectories decompose naturally into a fast 
component that contracts exponentially toward the slow manifold and a slow com- 
ponent which obeys the reduced system dynamics on the manifold. In this sense, the 
fast variables become slaved to the observables, and knowledge of the slow manifold 
and of the reduced dynamics on it suffices to determine the full long-term system 
dynamics. 

The identification and approximation of slow manifolds is usually achieved by 
employing a reduction method. We briefiy list a number of these: Intrinsic Low 
Dimensional Manifold (ILDM), Computational Singular Perturbation (CSP), Method 
of Invariant Manifold (MIM), Approximate Inertial Manifold approaches, and Fraser- 
Roussel iteration, and we refer the reader to [U [8] for a more extensive listing. 



1.1 A class of iterative algorithms based on the zero-derivative 
principle 

In [1], we developed a class of iterative algorithms to locate slow manifolds for systems 
of Ordinary Differential Equations (ODEs) of the form 

v' = q{u,v), V G R^^f, ^ ■ ' 

where Ns -|- Nf = N. We treated the variables u as the observables (that is, as 
parametrizing the slow manifold we are interested in), and we assumed that there 
exists an Ng— dimensional, attracting, invariant, slow manifold C, which is given lo- 
cally by the graph of a function v = v{u). However, we emphasize that we did not 
need explicit knowledge of which variables are fast and which are slow, only that the 
variables u suffice to parametrize C. 

To leading order, the location of a slow manifold C is obtained by setting v' = 0, 
i.e., by solving q{u,v) = for v. Of course, the manifold defined by this equation 
is in general not an invariant slow manifold under the fiow of the full system (11. ip . 
This is only approximately true, since higher-order derivatives with respect to the 
(fast) time t are, in general, large on it. If one requires that v" vanishes, then the 
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solutions with initial conditions at the points defined by this condition depend only 
on the slow time to one order higher, as v' also remains bounded in the vicinity 
of this manifold. Similarly, demanding that successively higher-order time deriva- 
tives vanish, we obtain manifolds where all time derivatives of lower order remain 
bounded. The solutions with these initial conditions depend only on the slow time 
to successively higher order and thus approximate, also to successively higher order, 
solutions on the slow manifold. In other words, demanding that time derivatives of 
successively higher order vanish, we filter out the fast dynamics of the solutions to 
successively higher orders. In this manner, the approximation of the slow manifold C 
is improved successively, as well. This idea may be traced back at least to the work of 
Kreiss [H [TTl [12] , who studied systems with rapid oscillations (asymptotically large 
frequencies) and introduced the bounded derivative principle to find approximations 
of slow manifolds as the sets of points at which the derivatives are bounded (not 
large). The requirement here that the derivatives with respect to the (fast) time t 
vanish is the analog for systems f 1 1.1 1) with asymptotically stable slow manifolds. A 
similar idea was introduced independently by Lorenz in [13j, where he used a simple 
functional iteration scheme to approximate the zero of the first derivative, then used 
the converged value of this scheme to initialize a similar scheme that approximates 
the zero of the second derivative, and so on until successive zeroes were found to be 
virtually identical. See also [3] and [6] for other works in which a similar condition is 
employed. 

The elements of the class of iterative algorithms introduced in |3] are indexed 
by m = 0, 1, . . .. The m— th algorithm is designed to locate, for any fixed value of 
the observable mq, an appropriate solution, v = Vm{uo), of the (m + 1)— st derivative 
condition 



Here, the time derivatives are evaluated along solutions of f ll.lj) . In general, since 
condition fll.2p constitutes a system of Nf nonlinear algebraic equations, the solution 
Vm{uo) cannot be computed explicitly. Also, the explicit form of (11. ip . and thus also 
an analytic formula for the (m + 1)— st time derivative in Eq. (II. 2p . may be unavail- 
able {e.g., in Equation- Free or legacy code applications). In this numerical 
approximation for it has to be used. The m-th algorithm in the class generates an ap- 
proximation t># of Vm{uo), rather than fm(uo) itself, using either an analytic formula 
for the time derivative or a finite difference approximation for it. In either case, the 
approximation t>* to Vm{uo) is determined through an explicit functional iteration 
scheme, which we now introduce. 

The m = algorithm is defined by the map Fq : R^f — > R'^f 




(1.2) 



Fo{v) = v + H 
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where i?, which we label as the iterative step size, is an arbitrary positive number 
whose magnitude we fix below for stability reasons. We initialize the iteration with 
some value v^^^ and generate the sequence 

1^(^+1) =^q(^;W)|^^ 1^2,...}. 

The functional iteration is terminated when ||t>('"+^) — v^'''^] < TOLq, for some r > 1 
and a prescribed tolerance TOLq. The output of this zeroth algorithm is the last 
member, v^, of the sequence {v'^'^^^^}. 

Next, the m — 1 algorithm is defined by the map Fi : R^* R^f, 

F,{v)=v-H' (^^^ M, 

initiahzed with some value v^^\ It generates the sequence 

= r = 1,2,...} 

and the functional iteration is terminated when — v^^^\\ < TOLi, for some 

r > 1 and for a prescribed tolerance TOLi. The output of this first algorithm is the 
last member, vf, of the sequence 

The algorithm with general m is defined by the map : HJ^' — > R'^f, 

Fm{v) = V - i-nr^' (^^) {no, v), (1.3) 

seeded with some value v^^\ It generates the sequence 

' .('-+i)=F^(^;M)|r = l,2,...}. 



Here also, one prescribes a tolerance TOL^ and terminates the iteration procedure 
when -v^''^\\ < TOL„, for some r > 1. The output of this m— th algorithm is 

the last member of the sequence {v^"^^^^}, denoted by v*. 

As we show in this article, not only is the point {uq, v*) of interest for each indi- 
vidual m because it approximates {uo,v{uq)), but the entire sequence {(-Uojf*)}™ is 
also of interest because it converges to {uq, v{uo)) with a suitably convergent sequence 
{TOL^}. Hence, the latter point can be approximated arbitrarily well by members 
of that sequence, and the class of algorithms may be used as an integrated sequence 
of algorithms in which the output of the m— th algorithm can be used to initialize 
the (m + 1)— st algorithm. Of course, other initializations are also possible, and we 
have carried out the analysis here in a manner that is independent of which choice 
one makes. 
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This class of iterative algorithms was applied in to three examples: the two- 
dimensional Michaelis-Menten mechanism for which the one-dimensional slow mani- 
fold can be computed analytically to arbitrary precision, a five-dimensional nonlinear 
system with an explicitly computable two-dimensional slow manifold, and a seven- 
dimensional hydrogen-oxygen system with quadratic nonlinearities for which the man- 
ifold is not known explicitly. In the context of these three examples, we found that, 
for all of the values of m that we worked with, the m-th algorithm converged at an 
exponential rate to a fixed point. Moreover, in the two examples where the slow 
manifold can be computed, we also found that, for each algorithm, this fixed point is 
very close to the actual point on the slow manifold. In addition to showing the m-th 
algorithm converged for each m that we worked with, we also showed that the class of 
algorithms may be used in the integrated manner stated above. The closeness of the 
approximation to {uq,v{uo)) improved as we increased the order m of the algorithm 
used. 

More recently, van Leemput et al. [16] employed the first (m = 0) algorithm 
in the class to initialize Lattice Boltzmann Models (LBM) from sets of macroscopic 
data in a way that eliminates the stiff dynamics triggered by a bad initialization. 
They showed that the algorithm they derived converges unconditionally to a fixed 
point close to a slow manifold, and they used the algorithm to couple a LBM to a 
react ion- diffusion equation along the interface with good results |l7j . 

Our motivation for introducing this class of iterative algorithms in [1] was two- 
fold. First, we wanted a method that can be implemented in the context of legacy 
codes. In other words, we wanted this reduction method to be implementable even 
when one has no explicit form for the components p and q of the vector field, but only a 
black-box integrator (timestepper). This feature renders the method "equation-free" 
[TU] and makes its implementation possible in these settings. Second, it was essential 
for us that they preserve the user-specified value of the observables, say m = Mq, at 
each iteration. In this way, the output of the algorithm is an approximation of the 
point {uq,v{uq)) on the manifold C of that same value Uq of the observables. Also, in 
this way, the 'lifting' step in projective integration of ^ is naturally facilitated. 

It is worth noting that one really only needs to require that the time derivatives 
are sufficiently small, although we work with the zero-derivative condition (11.21) for 
definiteness. 
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1.2 Iterative algorithms based on the zero-derivative princi- 
ple for explicit fast— slow systems 



A central assumption that we made in |1] is that we work with systems (11. II) for which 
there exists a smooth and invertible coordinate change 

z = z{w) with inverse w = wlz), (1-4) 

where w = {u, v) and z = {x, y), which puts the system f ll.ip into the exphcit fast-slow 
form 

x' = f{x,y,e), xgR^% ,^ 

We emphasize that, in general, we have no knowledge whatsoever of the transfor- 
mation that puts system fll.ip into an explicit fast-slow form. Here, / and g are 
smooth functions of their arguments, the manifold C is transformed smoothly, and 
det{Dyg)Q{z) = det{Dyg{z,0)) 7^ on the manifold £[o] = {z\g{z,0) = 0} (on which 
the dynamics reduce for e = 0), see also 

Due to the above assumption, it turns out to be natural to split the analysis of 
the accuracy and convergence of the functional iteration into two parts. In the first 
part, which we present in this article, we work directly on systems that are already 
in explicit fast-slow form (II. 5p . In the context of these systems, the accuracy and 
convergence analysis may be carried out completely in terms of the small parameter e. 
The system geometry - the slow manifold and the fast fibers transverse to C - makes 
the convergence analysis especially transparent. Then, in the second part, we work 
with the more general systems (11.10 . For these, the accuracy analysis proceeds along 
similar lines as that for this first part, with the same type of result as Theorem 12.11 
below. However, the convergence analysis is considerably more involved than that 
for explicit fast-slow systems. For these general systems, one must analyze a series 
of different scenarios depending on the relative orientations of (i) the tangent space 
to £, (ii) the tangent spaces to the fast fibers at their base points on £, and (iii) 
the hyperplane of the observables u. Moreover, all of the analysis must be carried 
out through the lens of the coordinate change (11.41) and its inverse, so that it is less 
transparent than it is in part one. Part two will be presented as a subsequent article. 

As applied specifically to explicit fast-slow systems (II. Sp . the m— th iterative 
algorithm (ll.3p is based on the (m + 1)— st derivative condition, 

(^)(-o.y) = o. (1,6) 

In particular, for each m and for any arbitrary, but fixed, value of the observable 
Xq G K, one makes an initial guess for h{xQ) and uses the m-th iterative algorithm 
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to approximate the appropriate zero of this (m + 1)— st derivative, where the end 
(converged) result of the iteration is the improved approximation of h{xQ). 

For each m = 0, 1, . . ., the m— th iterative algorithm is defined by the map 

Fm{y) = y - i-Hr+' ^-^j (xo, y), (1.7) 

where H is an arbitrary positive number whose magnitude is 0{e) for stability reasons. 
We seed with some value y^^^ and generate the sequence 

{y(^+i)^F™(yM)|r = l,2,...}. (1.8) 

Here also, one prescribes a tolerance TOL^ and terminates the iteration procedure 
when - < TOL^ for some r > 1. The output of this m— th algorithm is 

the last member of the sequence {y^^~^^^, denoted by y*. 



1.3 Statement of the main results 

In this article, we first examine the m-th iterative algorithm in which an analytical 
formula for the (m + 1)— st derivative is used, and we prove that it has a fixed point 
y = hm{xo), which is 0{e'^~^^) close to the corresponding point h{x()) on the invariant 
manifold C, for each m = 0, 1, . . .. See Theorem 12. II below. 

Second, we determine the conditions on {Dyg)^ under which the m-th iterative 
algorithm converges to this fixed point, again with an analytical formula for the 
(m + 1)— st derivative. In particular, for m = 0, the iteration converges for all 
systems (11. 5p for which {Dyg)o is uniformly Hurwitz on £[o] and provided that the 
iterative step size H is small enough. For each m > 1, convergence of the algorithm 
imposes more stringent conditions on H and on the spectrum of {Dyg)^. In particular, 
if a{{Dyg)Q) is contained in certain sets in the complex plane, which we identify 
completely, then the iteration converges for small enough values of the iterative step 
size H, see Theorem 13. 1[ These sets do not cover the entire half-plane, and thus 
complex eigenvalues can, in general, make the algorithm divergent. 

Third, we show explicitly how the Recursive Projection Method (RPM) of Shroff 
and Keller [15j stabilizes the functional iteration for each m > 1 in those regimes 
where the iteration is unstable. This stabilization result is useful for practical imple- 
mentation in the equation- free context; and, the RPM may also be used to accelerate 
convergence in those regimes in which the iterations converge slowly. Alternatively, 
the Newton-Krylov Generalized Minimal Residual Method (NK-GMRES [9]) may be 
used to achieve this stabilization. 
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Fourth, we analyze the influence of the tolerance, or stopping criterion, used 
to terminate the functional iteration. We show that, when the tolerance TOL^ for 
the m— th algorithm is set to 0{e"^~^^), the output also satisfies the asymptotic 
estimate \\y* - h{xo)\\ = C»(£™+^). 

Finally, we extend the accuracy and convergence analyses to the case where a 
forward difference approximation of the (m + 1)— st derivative is used in the iteration, 
instead of the analytical formula. As to the accuracy, we find that the m-th iterative 
algorithm also has a fixed point y = hjn{xo) which is close to h{xQ), so that 

the iteration in this case is as accurate asymptotically as the iteration with the analyt- 
ical formula. Then, as to the stability, we find that the m-th iterative algorithm with 
a forward difference approximation of the (m + 1)— st derivative converges uncondi- 
tionally for m = 0. Moreover, for m = 1, 2, . . ., the convergence is for a continuum of 
values of the iterative step size H and without further restrictions on {Dyg)o, other 
than that it is uniformly Hurwitz on C[o], see Theorem 16.11 These advantages stem 
from the use of a forward difference approximation, and we will show in a future work 
that the use of implicitly defined maps Fm yields similar advantages. 

Throughout this article, we shall refer to some basic facts about the Ng— dimensional, 
slow, invariant, and normally attracting manifold C. As stated above, C is the graph 
of a function h, 

C = {ix,y)\xeK,y = hix)}, (1.9) 
for some set K. Here, the function h : K H^'^ satisfies the invariance equation 

g{x, h{x),e) - eDh{x)f{x, h{x),6) = 0, (1.10) 

and it is 0{e) close to the critical manifold, which is the graph of ho{x), uniformly 
for X & K. 

It is insightful to recast this invariance equation in the form 



(— D/i(x), JnJ (^(a;, /i(x), e) = 0, where G = 




(1.11) 



which reveals a clear geometric interpretation. Since C corresponds to the zero level 
set of the function —h{x) + y hj Eq. (11.91) . the rows of the Nf x N gradient matrix 
{—Dh{x), Jnj) form a basis for NzC, the space normal to the slow manifold at the point 
z = {x, h{x)) G C Thus, Eq. (11. lip states that the vector field G is perpendicular to 
this space and hence contained in the space tangent to the slow manifold, TzC. 
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2 Existence of a fixed point hm{xo) and its proxim- 
ity to h{xo) 



We rewrite the map F^, given in Eq. fll.7p . as 

Fmiy) = y - Lmixo,y), (2.1) 

where the function '■ —>■ is given by 

LM^i-Hr^'l^^^y^), for any m = 0,1,... , (2.2) 

where z = {xo,y). The fixed points, y = /im(xo), of are determined by the 
equation 

Lm{Xo,hmiXo)) = 0, 

that is, by the (m + 1)— st derivative condition (11.61) . The desired results on the 
existence of the fixed point hm{xo) and on its proximity to h{xo) are then immediately 
at hand from the following theorem: 

Theorem 2.1 For each m = 0,1, . . ., the (m + l) — st derivative condition lil.6\) . 

L^(x, y) ^ i-nr^' j {x, y) = 0, (2.3) 

can be solved for y to yield an Ns— dimensional manifold Cm which is the graph of a 
function '■ K ^ R'"^' over x. Moreover, the asymptotic expansions of km and h 
agree up to and including terms ofO{e"^), 

m 

i=0 i=l 

This theorem guarantees that, for each xq G K, there exists an isolated fixed point 
y = hm{xo) of the functional iteration algorithm. Moreover, this fixed point varies 
smoothly with Xq, and the approximation (xq, hmixo)) of the point (xq, h{xo)) on the 
actual invariant slow manifold is valid up to 0{e"^^^). 

The remainder of this section is devoted to the proof of this theorem. We prove 
it for m = and m = 1 in Sections 12.11 and 12.21 respectively. Then, in Section 12. 3[ 
we use induction to prove the theorem for general m. 
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2.1 Proof of Theorem 12.11 for m = 



We show, for each x E K, that Lq{z) has a root y = hQ{x), that Hq hes 0{e) close to 
h[o]{x), the corresponding point on the critical manifold, and that the graph of the 
function ho over K forms a manifold. 

For 171 = 0, definition fl2.2p . the chain rule, and the ODEs (11.51) yield 

Lo = -Hy' = -e-'Hg. (2.4) 

Substituting the asymptotic expansion y = ho{x) = J2i=o^^^o,i{x) into this formula 
and combining it with the condition Lq = 0, we find that, to leading order, 

g{x, ho^o{x),0) = 0, 

where we have removed the 0{1), nonzero, scalar quantity —H/e. In comparison, the 
invariance equation (11.101) yields 

(7(x,/i[o](x),0) =0, (2.5) 

to leading order, see Eq. ( 1A.2I) in Appendix A. Thus /to,o can be chosen to be equal 
to /i[o], and Lo{z) has a root that is (9(e)— close to y = h{x). 

It remains to show that the graph of the function Hq is an Ng— dimensional 
manifold Cq. Using Eq. (12. 4p . we calculate 

{DyLo) = -e-'H{Dyg), 

where all quantities are evaluated at {x,ho{x),e). Moreover, 

(DyLo) (x, ho{x)) = -E~'H (Dyg)^ + 0(e), 

with (•)o = (■)(x, /io,o(a^), 0) = (■)(x, /i[o](a;), 0), since /iq.o = ^[o]- Thus, the Jacobian 
{DyLQ){x, ho{x)) is non-singular for < e ^ 1, because H = 0{e) by assumption 
and because det{Dyg)Q ^ 0, see the Introduction. Therefore, we have 

det (DyLo) {x, ho{x)) ^ 0, for all x E K, 

and hence Co is a manifold by the Implicit Function Theorem and [TH Theorem 1.13]. 
This completes the proof of the theorem for the case m = 0. 

2.2 The proof of Theorem 12.11 for m = 1 

In this section, we treat the m = 1 case. Technically speaking, one may proceed 
directly from the m = case to the induction step for general m. Nevertheless, we 
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find it useful to present a concrete instance and a preview of the general case, and 
hence we give a brief analysis of the m = 1 case here. 

We calculate 

Li = {-Hfy" = -H{-Hy'y = -HL', = ~e-^H{DM)G. 

Using the ODEs (ll.Sp and Eq. (12. 4p . we rewrite this as 

Li = {-s-^hY [e{D,g)f + {Dyg)g] . (2.6) 

We recall that the solution is denoted hy y = hi{x) and that we write its asymptotic 
expansion as hi{x) = J2i=o^^^hii^)- Substituting this expansion into Eq. (12. 6p and 
recalling that H = 0{e), we obtain at 0{1) 

L, = i-E-'HfiDyg)^go + 0{e), 

where (■)o = {■){x, hi^o{x),0). Hence, y = h[o-\{x) is a root of Li to leading order by 
Eq. (12. 5p and det^Dyg)^ ^ 0, and therefore hi^ can be selected to be equal to /i[o]. 

At 0{e), we obtain 

{-e-'Hy{Dyg), [{Dyg),\D,g)^h + (^,^?)o/^i,i + {D,g)^] = 0, (2.7) 

where we used the expansion 

g{x, h, e)=go + e [{Dyg)oh^i + {D,g)o] + 0{e^) 

and that go = g{x, hi Q, 0) = g{x, /i[o], 0). Differentiating both members of the identity 
g{x, /i[o](x), 0) = with respect to x, we obtain 

{D,g)o + {Dyg)o{Dh[o]) = 0, 

whence {Dyg)Q^{Dxg)Q = —Dh[Q\. Removing the invertible prefactor {—H/ey{Dyg)Q, 
we find that Eq. (12. 7p becomes 

-(D/i[o])/o + {Dyg)oh^i + {D,g)o = 0. 

This equation is identical to Eq. (1A.3P in Appendix A, and thus hi i = h[i]. Hence, 
we have shown that the asymptotic expansion of hi{x) agrees with that of h{x) up 
to and including terms of 0{e), as claimed for m = 1. 

Finally, the graph of the function hi forms an Ng— dimensional manifold Ci. This 
may be shown in a manner similar to that used above for Co in the case m = 0. This 
completes the proof for m = 1. 
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2.3 The induction step: the proof of Theorem [2JJ for general 

m 



In this section, we prove the induction step that estabhshes Theorem 12.11 for all m. 
We assume that the conclusion of Theorem 12.11 is true for m and show that it also 
holds for m + 1, i.e., that the condition 

[{D,L^){x,y)]G{x,y,e) = Q (2.8) 

can be solved for y to yield y = hm+i{x), where 

m+1 

To begin with, we recast the (m + 1)— st derivative condition Eq. (12.31) in a form that 
is reminiscent of the invariance equation, Eq. (II. lip . Let m > be arbitrary but 
fixed. It follows from definition (I22D, Eq. (|rTTil . and Eq. (O) that 

i„ . (i-Hr^) ^ -h'-^ ^ -e-^H(DMG. (2.9) 

Therefore, the (m + 1)— st derivative condition (12.31) can be rewritten in the desired 
form as 

(D,L,„_i) G = 0, (2.10) 
where we have removed the 0(1), nonzero, scalar quantity —H/e. 

The induction step will be now be established using a bootstrapping approach. 
First, we consider a modified version of Eq. (12.81) . namely the condition 

[(D,L„,) (a;, h^{x))] G(x, e) = 0, (2.11) 

in which the matrix DzLm is evaluated on Cm (already determined at the m— th 
iteration) instead of on the as- yet unknown Cm+i- This equation is easier to solve for 
the unknown ?/, since y appears only in G. We now show that the solution y = hm+i (x) 
of this condition approximates h up to and including 0{e"^~^^) terms. 

Lemma 2.1 The condition Eq. i[2.11\) can be solved for y to yield 

m+1 

l/ = /;,™+i(x) = ^e*/i[,](x) + C(£'"+^), forallxeK. (2.12) 

Then, with this first lemma in hand, we bootstrap up from the solution y = hm+i of 
this modified condition to find the solution y = hm+i of the full (m+ 1)— st derivative 
condition, Eq. (I2.10p . Specifically, we show that their asymptotic expansions agree 
up to and including terms of C(e'"~'"^), 



12 



Lemma 2.2 The condition i\2. 8\) . can be solved for y to yield 



m+l 



y = h^+i{x) = e"hm+iA^) + for all x e K. 



i=0 



Given these lemmata - the proofs of which are given in appendix [B] - Theorem 12.11 
follows directly. 



3 Stability analysis of the fixed point hm{xo) 

In this section, we analyze the stability type of the fixed point y = hm{xo) of the 
functional iteration scheme given by Fm{y)- To fix the notation, we let 

a{Dyg\ ={\, = Xe,R + % \i = = A,,k(1 + ^ tan^,) : £ = 1, . . . , Nf } (3.1) 

and remark that normal attractivity of the slow manifold implies that X^^r < 
(equivalently, n/2 < Oi < 37r/2) for all £ = 1, . . . , Nf. Then, we prove the following 
theorem: 

Theorem 3.1 For each m = 0, 1, . . the functional iteration scheme defined by 
is stable if and only if the following two conditions are satisfied for a// £ = 1, . . . , Nf .■ 



n o I I f2m + Ak + l 2m + Ak + 3 \ ^ \ f tt SttX 



k=0,...,m 

and 



(3.2) 



0<H < //^^^ = [2cos((m + l){ee - Tf))]'^^""^'^ . (3.3) 

In particular, if Xi, . . . , An^ are real, then the functional iteration is stable for all H 
satisfying 

H < i/"^^" = 2^/^™+^) „ ^ „ . (3.4) 

WDygh 

The graphs of the stability regions for m = 0, 1, 2, 3 are given in Figure [H 

We now prove this theorem. By definition, hm{xo) is exponentially attracting if 
and only if 

a((DFj(/i^(xo)))cB(0;l), (3.5) 



13 



where B(0; 1) denotes the open ball of radius one centered at the origin. To determine 
the spectrum of {DFm){hm{xo)), we use Eq. (12.11) and Lemma [B. II to obtain 



{DFJ {y) 



hi - {DyLm) {xo,y) 

/n, - {-e-'H{Dyg){xo,y,0)r''' + O (e, \\go{x,,y)\\) . 



Letting y = hm{xo) in this expression and observing that ||5'o(a;o) ^m(a;o))|| = ^'(e) by 
virtue of the estimate hm = hQ + 0{e) (see Theorem 12. ip and Eq. (12. 5p . we obtain to 
leading order 

{DFJ (hUxo)) = - {-e-^HDyg)'^^^ , (3.6) 

where = {xq, hm{xo)) and the notation (■)o signifies that the quantity in paren- 
theses is evaluated at the point (xq, h^{xo)) G £[o]- Finally, then, we find to leading 
order 



a {(DF^) (/i™(xo))) = {fie = l- (|A,| b-'H) 



l,...,Nf . (3.7) 



In view of Eq. (13. 7p . condition (13. 5p becomes 



l-{\Xe\e-'Hr^'e 



< 1, 



for all 



1 Nf. 



(3.^ 



Here, we note that higher order terms omitted from formula (13. 7p do not affect sta- 
bility for small enough values of e, because the stability region B(0; 1) is an open set. 
Next, we study the circumstances in which this stability condition is satisfied. This 
study naturally splits into the following two cases: 



Case 1: The eigenvalues Ai, . . . , An^ are real. This is the case, for example, when 
the fast part of system (11.50 corresponds to a spatial discretization of a self-adjoint 
operator. Here, 6e = tt for all i, and thus condition (13. 8p reduces to 

< {\Xi\e-^Hy'^^ < 2, for all£= l,...,Nf, 
which further yields Eq. (13. 4p . 



Case 2: Some of the eigenvalues Ai, 
Using Eq. (13.71) . we calculate 



Anj have nonzero imaginary parts. 



l + {\Xe\e-'Hr^' {\Xe\e-'H) 



. m+1 



2cos((m + l)(^£-7r)) 



This equation shows that is a convex quadratic function of Convexity 
implies that, if there exists a solution H^^^ > to the equation \fj,i\ = 1, then l/x^l < 1 
for alio < H < Hf''^^. Plainly, = 1 implies 



r 



2cos((m + l)(0^ - tt)) = 0, 



14 



which yields condition (13.31) . Further, the condition that if™'^^, . . . , H^^^ be real and 
positive translates into condition (13. 2p . This completes the proof of Theorem 13.11 

For later comparison to the results of numerical simulations, it is useful to write 
formula (13. 3p explicitly for the first several values of m. For m = 0, formula (13.31) 
becomes 

H^'^'' = -^2cos^^, 

see Figure m We note that H^^"" > for all 9i G (7r/2, 37r/2), and thus the fixed point 
ho{xo) is stable for all < < i/"^'^^, where if™''^ = min<.(iff ^^). 

For m = 1, formula (13. 3p becomes 

^r' = ^V2cos(2^,), 

see Figure [H We see that, on (7r/2, 37r/2), H^^^ > only if 9i lies in the subin- 
terval (37r/4, 57r/4). Therefore, the fixed point /ii(xo) is stable if and only if (i) 
Oi e (37r/4, 57r/4), for all £ = 1, ... , Nf, and (ii) 0< H < H"'^'' = min^(if7^^). 

For m = 2, formula (13. 3p becomes 

H^^"" = -^[2cos(3e^)]^/^ 

see Figure [H Here also, H^^^ > on (7r/2,37r/2) only if 6i lies in the subinterval 
(57r/6, 77r/6). Thus, /i2(a^o) is stable if and only if (i) 9i G (Svr/G, Tvr/G), for all 
£ = 1, . . . , Nf, and (ii) <H < /f™^^ = min^(i7f ^^). 

For m = 3, formula (13. 3p becomes 

^^' = ^[2cos(4^,)]l/^ 

see Figure [H We observe that, on (7r/2, 37r/2), H^"^^ > only if 6*^ lies in the 
subdomain (7r/2,57r/8) U (Tvr/S, 97r/8) U (IItt/S, 37r/2). Therefore, the fixed point 
h{xo) is stable if and only if (i) 6^ e {n/2, Stt/S) U (Ttt/S, Qtt/S) U (llvr/S, 37r/2), for 
all £ = 1, ... , Nf, and (ii) 0< H < i/™^^ = miniiH^^''). 



4 Stabilization of the algorithm using RPM 

In the previous section, we saw that, for any m > 1, the m— th algorithm in our 
class of algorithms may have a number of eigenvalues that either are unstable or 
have modulus only slightly less than one. In this section, we demonstrate how the 
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Jl/2 Jt 3Jt/2 Jt/2 3jt/4 n 5lt/4 3lt/2 

m=2 m=3 




Figure 1: if™"^ as a function of 6^ G (vr/2, 37r/2), for m = 0, 1, 2, 3. H'^^^ is measured 
in units oi e/ \ \(\. The eigenvalue /i£ is stable for all < < H^°-^ . 



Recursive Projection Method (RPM) of Shroff and Keller [15] may be used to stabilize 
the algorithm or to accelerate its convergence in all such cases. 

For the sake of clarity, we assume that {DFm){hm{xo)) has M eigenvalues, la- 
belled {/ii, . . . , fipi}, that lie outside the disk B(0; 1 — 5), for some small, user-specified 
6 > 0, and that the remaining Nf — M eigenvalues {/im+i, • • • , /^Nf } lie inside it. We 
let P denote the maximal invariant subspace of {DFm){hm{xo)) corresponding to 
{fii, . . . ,/UAf} and P denote the orthogonal projection operator from R^f onto that 
subspace. Additionally, we use Q to denote the orthogonal complement of P in R^f 
and Q = — P to denote the associated orthogonal projection operator. These 
definitions induce an orthogonal direct sum decomposition of R'"^'* , 

RNf = p Q = pRNf ^ gj^Nf^ 

and, as a result, each y G R'^* has a unique decomposition y = p + q, with p = Py G P 
and q = Qy G Q. The fixed point problem y = F^iy) may now be written as 

p = PF,n{p + q), (4.1) 
q = QF,n{p + q). (4.2) 

The fundamental idea of RPM is to use Newton iteration on Eq. (14.11) and func- 
tional iteration on Eq. fl4.2p . In particular, we decompose the point y'-^^ (which was 



16 



used to generate the sequence {y^^~^^^ in Eq. (11. 8p ) via 

y(i)=p(i) + g(i) = Py(i) + g^(i). 

Then, we apply Newton iteration on Eq. (14. ip (starting with p'^^^) and functional 
iteration on Eq. (14. 2 p (starting with g*-^-*), 

pir+i) = p(r) ^ ^j^^ _ pi^DFUP^^^ + g(^)))P] ' PF„(p('^) + gM), . . 

The iteration is terminated when Hy^'''''^-' — y*-^^!! < TOLm, for some r > 1, as was also 
the case with functional iteration. 

Application of Theorem 3.13 from [15] directly yields that the stabilized (or 
accelerated) iterative scheme (14. 3 p converges for all initial guesses y^^^ close enough 
to the fixed point hm{xo), as long as 

1 i a{P{DFm{hm{Xo)))P) = {flu fiM}. 

In our case, this condition is satisfied for all H > 0, because the fact that C is 
normally attracting implies that each eigenvalue Xe of DyQ is bounded away from 
zero uniformly over the domain K on which the slow manifold is defined. Thus, the 
iteration scheme (14.30 converges. 



5 Tuning of the tolerance 

In this section, we estabhsh that, for every m = 0, 1, . . ., ||?/* — /i(xo)|| = 0{e"^^^) 
whenever TOL^ = 0{e"^~^^). The value returned by the functional iteration is within 
the tolerance of the point on the true slow manifold for sufficiently small values of 
the tolerance. 

The brunt of the analysis needed to prove this principal result involves showing 
that, for these small tolerances, ?/* is within the tolerance of the fixed point, hm{xo)- 
The desired principal result is then immediately obtained by combining this result 
with the result of Theorem 12. H where it was shown that \\hm{xo) — h{xo)\\ = 0{e"^~^^). 

We begin by observing that 

\\y* - hmixo)\\ < \\y* - y^^'^W + \\y^''^ - h„,{xo)\\, for any r > 0, 

by the triangle inequality. The first term is 0{e^^^) by definition, as long as r is 
chosen large enough so that the stopping criterion, < TOL^, is satisfied. 

As to the second term, we may obtain the same type of estimate, as follows: First, 
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where we used Eq. (12.11) . and hence 

Second, is invertible in a neighborhood of its fixed point, by the Imphcit Function 
Theorem, because the Jacobian of Lm{xQ, ■) at hm^xo) is 

{DyLm) {zm) = {-s'^HDyg)'^^^ , 

by Eq. fl3.6l) . and det{Dyg) ^ since £[o] is normally attracting. Third, by combining 
these first two observations, we see that 

where denotes the local inverse of Lm{xo,-). Fourth, and finally, by expanding 
around zero, noting that L^{0) = hm{xo), and using the triangle inequality, we 
obtain 

\\y^'^ - h^ixo)\\ < ||(D,L-i)(0)|| llyW - + O (llyW - y('-+^)f ) . 

Recalling the stopping criterion, we have therefore obtained the desired bound on the 
second term, as well, 

WV^''^ - hUxo)\\ < II {DyL;^') (0)||TOL„ + O ((T0L„)2) . 

Hence, the analysis of this section is complete. 



6 The effects of differencing 



In a numerical setting, the time derivatives of y are approximated, at each iteration, 
by a differencing scheme, 

" ^ (^'"^'^^ ^-(^o,y) and H>0. 

In this section, we examine how the approximation and convergence results of Sec- 
tions [2H5] are affected by the use of differencing. We choose forward differencing, 

(A™+iy) (z) = J2{-ir^'-' (""J ) <P'i^;^H), (6.1) 

where (j){z] t) is a (numerically generated) solution with initial condition z, for con- 
creteness of exposition and where H is a. positive, 0{e) quantity. Also, forward 
differencing is directly implementable in an Equation-Free or legacy code setting. 
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By the Mean Value Theorem, 



H 



m+l 

1 

7] 



m+l 



z - 



m + l 
2 

m + 1 ^ 
2r] 



H 



-'m+l 



m+2 



d'^+^y 
~dF^ 



{<P{z-t)) 



(6.2) 



where rj = H/H > is an 0{1) parameter available for tuning and (f){z; t) is the point 
on the solution (f){z; t) at some time t G [0, (m + !)-&]. Thus, for the m— th algorithm, 
the approximation of d'^'^^y / df^^^ by the above scheme corresponds to generating 
the sequence {|/^''^|r = 1, 2, . . .} using the map 



Fm{y) = y - Lm{z), z = {xo,y), 



where 



Lm{z) = {-vr+' (A-+iy) (z) = LUz) 
Therefore, by Eq. (16.21) . 



m + l 

27] 



-'m+l 



{<P{z;i)). 



(6.3) 
(6.4) 



Fm{y) = Fm{y) + Lm+l{(t){z\t)). 

2ri 



Remark. For convenience in the analysis in this section, we take the flow (p to be 
the exact flow corresponding to Eq. (11. 5p . The analysis extends directly to many 
problems for which only a numerical approximation of (j) is known. For example, if 
the discretization procedure admits a smooth error expansion (such as exists often 
for flxed step-size integrators in legacy codes or in the Equation- Free context), then 
the leading order results still hold, and the map obtained numerically is sufficiently 
accurate so that the remainder estimates below hold. In particular, given a p-th order 
scheme and an integration step size h, it suffices to take h = 0{e) to guarantee that 
the error made in using the numerically-obtained map is 0{eP). Of course, with 
other integrators, one could alternatively require that the timestepper be 0{e"^~^'^) 
accurate, i.e., of one- higher order of accuracy. 

6.1 Existence of a fixed point hm{xQ) of the map Fm 

In this section, we establish that the map Fm has an isolated flxed point y = hm{x) 
which differs from hm{xQ) (and thus also from h{xo), by virtue of Theorem 12.11) only 
by terms of 0(e™+i). 

The flxed point condition Fm{xo,y) = y may be rewritten as 

<jji -|- 2 

= Lm{xo, y) = Lm{xo, y) — Lm+i(0(xo, y; t)), (6.5) 
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where we combined Eqs. (16.31) and (16.41) . In order to show that Fm has an isolated 
fixed point hm{xQ) which is 0(£™+^)— close to hm{xo), we need to establish the validity 
of the following two conditions. 



(i) The second term in the right member of Eq. (16. 5p satisfies the asymptotic esti- 
mate 

||L„+i(0(2„;t))|| = 0(e™+^), where = {xo,hmixo)). (6.6) 



(ii) The Jacobian of satisfies 



det [DyLm] (zm) 7^ and [DyLm] (zm) ^ = 0{1). 



(6.7) 



Let us begin by examining the term Lm+i{(f){zm', i))- Let (x, y) = (j){zm', i). Then, 
we may write 

Lm+ii(pizrr,;i)) = (x, y) - L^+i (£ , (x) ) , 

because Lm+i{-, hm+i{-)) = by the definition of L^+i and hm+i- Hence, 

\\Lm+i{4>{zm;i))\\ < \\{DyLm+i){x,hm+i{x))\\ \\y - hrn+i{x)\\ + O {\\y - hm+i{x)\\'^) . 

(6.8) 

Now, \\{DyLm+i){x, hm+i{x))\\ is 0{1) by Lemma iB.ll Next, the triangle inequality 
yields 

\\y - hm+i{x)\\ < \\y - h{x)\\ + \\h{x) - h.m+i{x)\\. 

The first term in the right member remains (9(e™'+^) for all times i G [0, (m + 1)H)]. 
Indeed, the initial condition is 0{e"^~^^)-c\ose to the normally attracting manifold 
C. Thus, the Fenichel normal form [7] guarantees that the orbit generated by this 
initial condition remains 0{e"^~^^)-close to C for 0{1) time intervals. The second 
term in the right member is also (^(e™^^), by Theorem 12.11 Thus, \\y — hm+i{x)\\ 
is also 0{e"^~^^). Substituting these estimations into inequality (16. Sp . we obtain that 
||Lm+i(0(2;m; ^))|| is 0{e"^~^^) and condition (16. 6p is satisfied. 

Next, we determine the spectrum of {DyLm){zm) to leading order to check con- 
dition (16. 7p . We will work with the definition of A™"*"^?/, Eq. (16.10 . rather than with 
formula (16. 2p which involves the unknown time t. Combining Eqs. (16.10 and (16. 3p . 
we obtain 

m+1 



7] 



m+1 



E 

^=0 



m + 1 



{-i)%y{z-m). 
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Differentiating both members of tliis equation with respect to y, we obtain 

m+l , I 1 \ 

{d,L^ {z)=r'~'Y.y I ) i-^YiDy<py){z-AH). (6.9) 

Next, {Dy(j)y){zm;t) = Q'^''^l'^)^^y9)o ^ leading order and for all t of 0{e) by standard 
results. Since IH = 0{e) for all £ = 0, 1, . . . , (m + 1), we may use this formula to 
rewrite Eq. (16. 9p to leading order as 

/ 1 \ f 4-1 

n—n \ / 



i=0 

Hence, 

a i i DyL^ ) {zj ] = I r?"+i ( 1 - e^^^/^ ^ 



l,...,NfL (6.10) 



where Zm = [xq, hm{xo)). This leading order formula for the elements of the spectrum 
shows that {DyLm){zm) is 0{1) and non-degenerate for all positive 0{e) values of H 
and H. Thus, condition (16.71) is also satisfied. 



6.2 Stability of the fixed point hm{xo) for 77 = 1 

In this section, we determine the stability of the fixed point hm{xo) under functional 
iteration using Fm in the case that H = H. Our results for H = H are summarized 
in the following theorem. The general case H ^ H is treated in the next section, and 
the main result there is given in Theorem 16. 2[ 

Theorem 6.1 Fix rj = 1. The functional iteration scheme defined by Fq is uncon- 
ditionally stable. For each m = 1,2, . . the functional iteration scheme defined by 
Fm is stable if and only if, for each i = 1, . . . , Nf, the pair {H, 6^) lies in the stability 
region the boundary of which is given by the implicit equation 

j=i fc=i ^ ^ ^ 

m+l , I 1 \ 2 

+ E( ) '''''''''' H, = -\RH/e>Q. (6.11) 

fc=i ^ ^ 

Here, the branch 0/ arctan is chosen so that 6e G (7r/2, 37r/2). In particular, if 
Ai, . . . , are real, then the functional iteration is unconditionally stable. If at least 
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one of the eigenvalues has a nonzero imaginary part, then a sufficient and uniform 
(in $1, ... , 6^ J condition for stability is that 

H > where Hs(l) = -In (2^/(™+^) - l) > 0. (6.12) 

The stability regions for various values of m are plotted in Figure [31 

Following the procedure used in Section [3l we determine a{{DFm){hm{xQ))) and 
examine the circumstances in which the stability condition 



(6.13) 



ai^i^DF^J {hUxo))) CB(0;1) 
is satisfied. Equation (16. 3p yields 

and thus also 

{fie} = cr ( (DyFm) (h^{xo))) = l- a ( [DyLm] (xq, h^{xo) 



Since hm{xo) differs from hm{xo) only at terms of 0{e'^^'^), {DyLm){xo,hm{xo)) also 
differs from {DyLm){xo, hm{xo)) only at terms of C(e™"''^). Thus, Eq. ( ]6.10p yields, 
to leading order and for £ = 1, . . . , Nf, 



= 1 _ (1 _ ^^.H/s^"^^^ = J2( k j i-lf^'e''^''^'. (6.14) 

k=i ^ ^ 

Recalling Eq. (13.11) and defining = —Xi^nH/e, we rewrite Eq. (16.141) in the form 



m+l ^ 
k=l ^ 



k 



(6.15) 



The stability condition (16.131) becomes, then, 

m+l 



\M 



E 

k=l 



m+l 
k 



< 1, for alH= l,...,Nf. (6.16) 



As in Section [3], we distinguish two cases. 
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Case 1: All of the eigenvalues of {Dyg)o are real. Then, Oi = tt for all £ = 

1, . . . , Nf, and hence Eq. (16.151) becomes 




Thus, the spectrum of {DyFrn){hm{xo)) is contained in (0,1) for all positive 0{e) 
values of H. Equivalently, the fixed point hm{xo) is unconditionally stable for these 
values of H. 

These results may be interpreted both in the context of the m-th iterative al- 
gorithm for each fixed m, as well as in the context of using the algorithms as an 
integrated class. In particular, for each fixed m, the rate of convergence to the fixed 
point of the m-th algorithm increases as H increases. Also, for any fixed iterative 
step size H, the rate of convergence of the m-th algorithm to its fixed point decreases 
as the order, m, of the iterative algorithm increases. This information is important 
for determining how large an m one should use, especially when using the algorithms 
as an integrated class. 



Case 2: Some of the eigenvalues of {Dyg)o have nonzero imaginary parts. 

When this is the case, some of the eigenvalues may be unstable for certain values of 
H. Figure [2] demonstrates this: in it, we have drawn the complex eigenvalue fii for 
various values of H and for m = 0, 1, 2, 3. Plainly, fi£ is unstable for m > and for H 
small enough, as > 1. We determine the stability regions in the (6'^, if^)— plane 
as functions of m. 

First, we derive the uniform bound (I6.12p . Using formula (16.151) . we calculate 

\M<J2i )e-'^^ = (l + e-^0"^+'-l, (6.17) 

fc=i ^ ^ 

and thus \fii\ < 1, for all > Hs{l). Recalling that Hi = —Xi^jiH/e, we conclude 
that all of the eigenvalues fii lie in the unit disk (equivalently, the m— th algorithm 
is stable) for all 0{e) values of H greater than eHs{l)/ min^ \ \i^r\, irrespective of the 
values of ^1 , ... , 6^^ . This is demonstrated in Figure [31 

Next, we derive formulae which describe exactly the stability regions. For m = 0, 
Eq. (I6.12P yields -f^s(l) = 0. Thus, < 1 for all positive 0{e) values of H and 
for all £ = 1, ... , Nf. As a result, the fixed point /io(a^o) is unconditionally stable for 
positive, 0{e) values of if, see also Figure [31 

For m = 1, Eq. (16.151) becomes 

p. _ r^ -Hi{l+itaiiei) _ -2Hi{l+itaiiei) 
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m=0 



m=1 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 

m=2 m=3 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 



Figure 2: The eigenvalue fii for values of H between zero and lOOe. The thick line 
denotes the boundary of the stability region {i.e., the unit circle). The eigenvalue Xi 
was taken to be — 1 + i for each one of the graphs. The arrow points to increasing 
values of H. 



Writing fi£ for the complex conjugate of fii, then, we calculate 



fie fie = 4e"^^* - 4e"^^* cos{H(, tan6'<.) + e" 



■AHf 



(6.18) 



Using this formula, we recast the stability condition fl6.16p into the form 



4e- 



-2Hf 



4e"^^'^ cos(if£tan6'£) + e 



< 1. 



In particular, the boundary of the stability region can be obtained by equating the 
expression in the left member of this inequality to one and solving for 6e, to obtain 



Of = arctan ( H, ^ 



arccos 



1 

-e 



+ e 



1 



He 

4" ' ' 4 



-e 



+ 2k7r 



Here, k E Z and the branch of arctan is chosen so that Be G (7r/2, 37r/2). We have 
plotted the stability region in Figure [31 We also note here that the boundary of the 
stability region close to 7r/2 and to 37r/2 has fine structure, see Figure HI 

For a general value of m, the stability condition fl6.16p is 

m+l 



E 

k=l 



m+l 
k 



< 1, for all £ = 1, 



Nf. 
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ID 




Figure 3: The regions of H for which < 1 as functions of 6e G (vr/2, 37r/2). White 
corresponds to stabihty {\fii\ < 1) and black to instabihty > 1). H is measured 
in units of e/ \Xe^R\- The angle 6e takes values on (7r/2, 371/2) and the black horizontal 
line corresponds to the uniform bound -^^(l) of Eq. (16. 121) . 



m + 1 
J 

m + 1 
j 



Now, using Eq. fl6.15p . we calculate 

m+l m+1 

= EE 

j=i k=i 

m+l j—1 

= 2EE 

3=1 k=l 

""^^ ^ m + 1 y ^_2fc//^ 
A; 



m + l 

m + l 
A; 



^ _ ]^ y'+fcg-(i+fe)-f^£' g«0 -fc)-^* tan 6 



.l)i+fce-0+fc)^^ cos ((j - k)Hetcinee) 



+E 

fc=i 



Equation (16. lip now follows directly. 
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0.503X 0.5097[ 0.51 67U 0.522z 0.528ir 0.5357t 



Figure 4: The fine structure of tlie stability region depicted in Figure [3] (witli m = 1) 
close to 7r/2. The exterior of the lobes is part of the stability region. 



6.3 Stability of the fixed point hm{xo) for rj 1 

In this section, we determine the stability of the fixed point hm{xo) for H ^ H. We 
define the function 

Our results are summarized in the following theorem. 

Theorem 6.2 Fix rj > 0. For each m = 0, 1, 2, . . the functional iteration scheme 
defined by F^ is stable if and only if, for each i = 1, . . . , Nf , the pair {H, Oi) lies in 
the stability region the boundary of which is given by the implicit equation 



j=l A; = l ^ / V / 

+2^-+i (^m+i - 1) 5^ ( + M (-l)^e-'=^*cos (km^nO, 



fc=l 

m+1 . , 1 \ 2 

+^2(m+i) ^ ( + M e-^^^^ + (ry-+i - 1)\ (6.20) 



26 



where He = —Xi^rH/e > 0. Here, the branch of arctan is chosen so that 6^ G 
(7r/2, 37r/2) . In particular: 

(i) Assume that \m{\e) = 0, for all £ = 1, . . . If < r] < 2^/^'^+'^\ then the 
functional iteration is unconditionally stable. If rj > 2^^^"^~^^\ then the functional 
iteration is stable if and only if 

(ii) Assume that at least one o/Im(Ai), . . . ,Im(ANf) is nonzero. If < t] < 2^/^"^~^^\ 
then a sufficient and uniform (in 9i, . . . , 9^ J condition for stability is 

H > JEh^^. (6.22) 

If rj > 2-'^/('"+-'^\ the functional iteration is unstable for any 9i, . . . , 6*^^ and for all 

H > (6.23) 

These results are demonstrated in Figures [5] and [61 

As in Section 16. 2[ we determine when the stabihty condition fl6.13p holds. The 
analogue of Eqs. (16.141) and (16.151) in this case is, to leading order and for £ = 1, . . . , Nf , 

jj^ = 1 _ 1 _ e^^-f^/M ^i_^m+i /^_^-H,(i+it^ne,)\ .(6.24) 

The stability condition (16.131) becomes, then, 

1 _ e-H,(l+*tan9,) j ^ g^jj £ = 1, . . . , Nf. (6.25) 

Here also, we distinguish two cases. 



Case 1: All of the eigenvalues of {Dyg)o are real. Then, 9e = n for all i = 

1, . . . , Nf, and hence Eq. (16.251) becomes 

/i^ = 1 -^'"+1(1 -6"^^)'"+^ 

Plainly, the condition fii < 1 is satisfied for all positive He and rj. Next, solving this 
equation for ri, we obtain an equation for the level curve fie = constant, 

1 _ e-H, 
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Figure 5: The stability region in the [r], Hfj— plane together with the level curves 
fii{ri, Hi) = —1 (thick curve), //^(r/, Hi) = (sohd curve in the middle), fii{ri, Hi) = 1 
(union of the two semiaxes). The dashed level curves to the right and left of the 
level curve fii = correspond to representative positive and negative values of fii, 
respectively. The eigenvalue fii is stable for all pairs {t], Hi) to the left of the level 
curve fii = —1. 



For < 1] < 2^/(™'+^) and for all 0{e) and positive values of H, we obtain fii > —1 (and 
thus the eigenvalue fii is stable), see Fig. [5l Therefore, a{{DyFm,){hm{xo))) C (—1, 1), 
and the fixed point hm{xo) is unconditionally stable. 

For 1] > 2^/(™+^), we obtain the condition < Hi < Hmiv), and Eq. 06.211) 
follows directly. Finally, we note that, for a fixed value of t] and as H oo, the 
spectrum clusters around 1 — ry™+^. Thus, the choice r/ = 1 is optimal in the sense 
that large values of H bring the spectrum closer to zero. 



Case 2: Some of the eigenvalues of {Dyg)^ have nonzero imaginary parts. 

In this case, some of the eigenvalues may become unstable for certain combinations 
of 7] and -ff, as our analysis in Section 16.21 also showed. 

First, we consider the case < < 2^/^'"+^) and derive the uniform bound (I6.22p . 
Using formula fl6.24p and working as in Eq. (16.171) . we estimate 



m+l 



:i+e 



-Hi\m+l 
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Figure 6: The stability regions in the (77, ^^)— plane for m = (left panel) and 
m = 1, 2, . . . (right panel). The eigenvalue fig is stable in region I, unstable in region 
II, and its stability type is dependent in region III. 



Hence 



< 



1+V 



m+l 



:i + e- 



-Hi\m+1 



for < ?7 < 1, 
for f] > 1. 



Combining these inequalities with the stability condition \fi£\ < 1, we obtain the 
sufficient condition Hi > Hm{ri), where Hm{v) is the uniform bound fl6.19p (see 
also Fig. [6]). Recalling that Hg = —Xi^rH/e, we conclude that, if condition (16.221) 
is satisfied, then a{{DyFm){hm{xo))) C B(0; 1), and hence the m— th algorithm is 
stable. 

Next, we consider the case i] > 2^/*^™+^) and derive the uniform bound (I6.23p . 
Equation (I6.24p yields 



|1 - P'el > V 



m+l 



^-Hi^iHi ta.nl 



m+l 



> 7]""+^ i 1 - e-"' 



m+l 



Thus, \l- ile\ > 2, for r] > 2^/^"+^) and He > H^nijl), and therefore 

> ||1 - A^l - 1| > 1, 

Hence, fie, is unstable. 
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Ti=1.4 




3it/2 



ri=1 .41 421 35623730950 



3JI/2 





1 



37t/2 



ri=1 .41 421 35623730951 



3JI/2 




3JI/2 



Figure 7: The stability region in the {rj, ^f^)— plane for m = 1 and for various values 
of rj. The last two values for rj are just below and just above the value 2^/^'""'"^) = a/2. 
White denotes stability and black denotes instability. 



Remark. Conditions (16.221) and fl6.23p may be interpreted by means of the fact 
that a{{DyF^)ihm{xQ))) clusters around 1 - r/™+i as ^ oo. For < < 2^'^'^+^\ 
there holds that — 1 < 1 — 7^"^+^ < 1. Thus, for H large enough, the eigenvalues are 
contained in the unit disk. On the contrary, 1 — r/'""'"^ < —1 for rf > 2^/^'""'"^), and 
thus the eigenvalues lie outside the unit disk for H large enough. 

Finally, formula (16.201) describing the stability region may be derived in a manner 
entirely analogous to that used to derive Eq. (16. lip . 



7 Conclusions and Discussion 



In this article, we characterized the accuracy and convergence properties of the class 
of iterative algorithms introduced in for explicit fast-slow systems (11.50 . The 
m-th member of the class corresponds to a functional iteration scheme to solve the 
{m + 1)— st derivative condition (II. 6p . We showed that this condition has an isolated 
solution, which corresponds to a fixed point of this m-th member and which is accurate 
up to and including terms of 0{e"^), see Theorem 12. 1[ Also, we derived explicit 
formulae for the domain of convergence of the functional iteration, both in the case 
where analytical formulae for the (m + 1)— st derivative are used (see Theorem 13.11) 
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and in the case where the (m + 1)— st derivatives are estimated through a forward 
difference scheme (see Theorem 16. ip . These convergence results are illustrated in 
Figures [T], [HI and HI Further, we demonstrated how the Recursive Projection Method 
may be used to stabilize the functional iteration in all cases when it is unstable or to 
accelerate its convergence in those cases where the convergence is slow. 

An extension of the analysis presented here to more general multiscale systems 
(11. ip will be presented in a subsequent article. The analysis of the accuracy of the 
(m + 1)— st derivative condition presented in Section [2] carries through, essentially 
(modulo a number of technicalities), in the more general well. The analysis of 

the stability of the functional iteration, on the other hand, is far more involved. The 
reason for that is that, although the hyperplane u = uq and the space tangent to the 
fast fibration over the slow manifold coincide to leading order for explicit fast-slow 
systems (II. 5p . this is not the case for the more general systems (II. ip . The absence of 
this feature makes the stability question for the functional iteration far more difficult 
to answer in the general case. 

In addition, we are in the process of generalizing the results of this article to other 
maps that may be used in the context of the functional iteration scheme developed 
in [3|. In particular, it is of interest to use maps which are implicitly defined (as 
opposed to the explicitly defined ones presented in [3] and in this article). Preliminary 
analytical results for m = and m = 1 indicate that one may construct functional 
iteration schemes based on implicit maps which not only retain the accuracy of the 
functional iteration scheme presented in this article but which are also unconditionally 
stable. Moreover, we think that this analysis may be extended to higher values of m, 
and we note that it is also possible to carry out the functional iteration with implicitly 
defined maps even when one only has a legacy code as a timestepper. 



A The one-higher-order proposition 

In this appendix, we state and prove a technical proposition - called the one-higher- 
order proposition - about the asymptotic accuracy of approximations of C given an 
approximation of the normal space to C This result is instrumental in the proof of 
the technical lemmas contained in the next appendix. 

We begin by recalling the useful formulation, Eq. (II. lip , of the invariance equa- 
tion that defines the function h{x), whose graph is the invariant, slow manifold C 
This formulation revealed that the matrix {—Dh{x),lN^) forms a basis for N^C, the 
space normal to the slow manifold at the point z = {x, h{x)) G C. 
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The function h{x) admits an asymptotic expansion in e, 

M-) = $^e'%](-), (A.l) 



where the coefficients i = 0,1, . . . , are determined by expanding asymptotically 
the left member of Eq. (ll.lOp and setting the coefficient of equal to zero to obtain 



9i - ^iDh[i]) fi^i^e = 0, ii = 0, 1, 



where the sum is understood to be empty for i = 0. The ffist few equations are 

^0 = 0, (A.2) 

{Dyg)oh[,] + {D,g)o - {Dh[o])fo = 0. (A.3) 

Here, Eq. ( ]A.2p is satisfied identically, Eq. ( 1A.3I) yields the coefficient h[i], and so on. 

The one-higher-order proposition, which we now state and prove, establishes a 
connection between the order in e to which a set N of row vectors approximates N^C 
and the order to which the solution ri{x) to the condition N G = approximates h. 



Proposition A.l Let N{x,e) be an Nf x N matrix with the property that its rows 
span N^i2 up to and including terms ofOle"^), for some m = 0,1, . . . . That is, N{-,e) 
is of the form 

N{;e) = c(-f2^'Dh[^i-)- E ^'M-),InX (A.4) 

V t=0 i>m+l / 

where C is a non-singular Nf x Nf matrix and Ri ^ Dh^, for i = m + l,m + 2, . . . , 
in general. Then, the condition 

N{x,e)G{x,y,£) = (A.5) 

can be solved for y to yield a function y = ri{x), the asymptotic expansion of which 
agrees with that of h{x) up to and including terms of 0{e"^'^^), 

m+l 

r^(x) = Y.^\{x) = J2e'h[^{x) + (A.6) 

i=0 i=0 



This proposition is called the one-higher-order proposition, because it states that the 
order to which ri{x) approximates the full slow manifold is of one higher than that to 
which N approximates the normal space. 
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Proof of Proposition lA.ll We recall that h{-) = Ei=o£^*^[i](-)5 by Eq. (lA.ip . 
and that h^q is determined from the ^(e*) terms of the invariance equation (II. lip . 
Similarly, rji is determined from the terms of Eq. (1A.5I1 . Thus, to establish 

Eq. (1A.6I) . it suffices to compare the terms of these two equations from 0{1) up 
through and including 0{e"^~^^) and to show that they are equal. 

First, for each i = 0, 1, . . . , m, the invariance equation (11.111) at 0{e'^) is 

i 

(-D/i[o],/n,) G, + J2 {-Dhie],0) G,^e = 0. (A.7) 

Second, to derive the 0{e^) terms for the condition NG = 0, we substitute the 
hypothesis flA.4p in Eq. flA.Sp and left-multiply by to obtain 

C-^NG= + C(£'"+^),/Nfj ^ = 0. (A.8) 

Thus, for each i = 0,1, . . . ,m, this condition at 0{e^) is 

i 

(-D/i[o],/n,) Gi + Y, {-Dhie],0) Gi^, = 0. 



1=1 



Plainly, this equation is identical to Eq. (lA.7p . Thus, r]i = h[i], for i = 0, 1, . . . 



m. 



Finally, we look at the 0{e''^~^^) terms of the two equations. Eq. flA.7p with 
z = m + 1 is 



{-Dhio],I^,) Gm+i + Yl {-Dh[i], 0) Gm+i-e + {-Dh[m+i],0) Gq = 0. (A.9) 
Also, Eq. dEH) at 0{e"'+^) is 

m 

(-D/i[o], In,) Gm+i + Y (-^^M' 0) Gm+i^e + {Rm+i, 0) Go = 0. (A. 10) 

e=i 

We note that Rm+i 7^ —Dh[m+i], in general. However, Gq = 0, since the terms ap- 
pearing in Eqs. (lA.9p - (lA.10p are evaluated at (x, r^o, 0) = (x, h[o], 0). Thus, Eqs. (1A.9P 
and (lA.lOp also agree, and hence rj^^i = h[m+i]- This completes the proof of the 
proposition. | 
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B Proofs of Lemmata 12.11 and 12.2 



In this appendix, we prove lemmata 12.11 and 12.21 characterizing the asymptotic accu- 
racy of the approximation to C obtained from the (m + 1)— st derivative condition 



Proof of Lemma 12. IL We write for (x, hm{x)) and z for {x, h{x)). The strategy 
is as follows: We will show that the rows of {DzLm){zm,s) span N^il up to and 
including terms of 0{e'^). Then, we will apply Proposition lA.ll to establish Eq. fl2.12p . 

The manifold Cm is the graph of the function h^i and thus it coincides exactly 
with the zero level set of the function —hm{x) + y. As a result, the rows of the Nf x N 
gradient matrix {—Dhm{x), form a basis for Nz^Cm- Second, the function hm{-) 
is defined through the (m+ 1)— st derivative condition Lm{-, hm{-),£) = 0. Therefore, 
Cm also coincides with (a connected component of) the zero level set of the function 
Lm{z, e). Thus, the rows of the Nf x N gradient matrix {DzLm){x, h„i{x),e) also form 
a basis for N^^/^m. It follows from the existence of these two bases that there exists 
a non-singular Nf x Nf matrix C such that 

{DzLm) {■,hm{-),e) = C{~Dhm{-)jN,) • (B.l) 



Next, the induction hypothesis implies that the asymptotic expansions of hm and 
h agree up to and including terms of 0{e"^), 



hm{-) = J2e'h[^{-)+0{e'^+'). (B.2) 



1=0 



Since the vector field is assumed to be sufficiently smooth, we may differentiate both 
sides of this equation with respect to x to obtain 

m 

Dhm{-) = Y,e'Dh^^{-)+0{e"'^^). (B.3) 

j=0 



Combining Eqs. (]B.1|) and (IB. 31) . then, we find 



V i=0 

This equation shows that the rows of {DzLm){x, hm{x),e) span N^£ up to and includ- 
ing terms of 0{e"^). Hence, application of the one-higher-order proposition. Proposi- 
tion A.l, completes the proof of this lemma. I 
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Before we proceed with the proof of Lemma 12.21 we prove the following result 
which will be needed therein. 



Lemma B.l For m = 0,1, . . ., for H = 0{e), and for a general point z = {x, y), the 
function is written as 

where the notation "{■)q{z)" stands for {■){z,0). The Jacobian DyL^ is written as 
{DyLm){z) = {-e-'H{Dyg)^)"'^' + 0{e,\\go{z)\\). (B.4) 

Proof. For this proof, we write (■)o instead of {■)o{z) for the sake of brevity. The 
proof is by induction on m. For m = 0, we recall Eq. (12. 4p . 

Lq = -e'^Hg, 

and hence, expanding g in powers of e, we find 

Lo = -e-^Hg^ + 0{e). 

This is the desired formula for Lq. Differentiating both members of this formula with 
respect to y, we obtain 

DyU = -e-^H{Dyg)^ + 0{e). 

This is the desired formula for DyL^. 

Next, we carry out the induction step for general m, namely we assume that 

Lm = {-e-^H)'^^\Dyg)^g, + 0{e,\\gor), (B.5) 

DyL^ = {~6-'HiDyg)^)"'^' + 0{e,\\goiz)\\) (B.6) 

and show that 

W = {-s-'H)"'^\Dyg)^'-'go + 0{E,\\gor). (B.7) 

DyL^+, = {-6-'H{Dyg)^)"'^' + Oie,\\goiz)\\). (B.8) 

By Eq. (Q, 

Lm+i = -e-'H{D,L^)G = -e-'H [e{a,L^,)f + {DyLrn)g] , 

Then, we substitute the induction hypothesis (IB.SP into this expression. Application 
of the differential operator {—H/£)[£{D.ji:-)f + {Dy-)g\ on the 0{e, Wg^W^) remainder 
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does not alter its asymptotic magnitude. Moreover, the term e{DxLm)f is 0{e) 
and, hence, can be absorbed also in the remainder. Therefore, we are left with the 
term {—H/e){DyLm)g- Substituting DyLm into this expression from the induction 
hypothesis (IB. 61) . we arrive at the desired formula (IB. 71) . 

Finally, we prove the leading order formula ( IB. 81) . First, we differentiate both 
members of the leading order formula (1B.7P with respect to y and use the product 



rule derivative to evaluate the right member. The second term from the product rule 
is precisely the leading order term in Eq. (IB.4p . The other term from the product 
rule, 



m 



may be absorbed in the remainder since it is linear in g^. Thus, we have obtained the 
desired formula (IB.SP and completed the proof of the lemma. I 



Proof of Lemma l2.2i We first use [21 Theorem 3] to establish that condition (12.101) 
has a solution y = hm+i{x) which is close to hm+i- According to that 

theorem, it suffices to show that 

((D,L„)(a;,/i„+i(x),£)) G{x,hrn+i{x),e) = 0{e^+'). 

By the definition of hm+i, 

iiD,Lm)i; hmi-),e)) Gi; hm+i{-),e) = 0. 

Thus, we may write 

{D^Lm){-, hm+i{-),e)^ G{-, hm+i{-),£) 

{D,Lm){-~hm+i{-),e) - {D,Lm){-,hm{-),e)\ G{-~hm+i{-),e). (B.9) 

Next, we have the following estimates of the asymptotic magnitudes of the two 
terms in the right member of Eq. (IB. 91) : 

m+l 

by Lemma 12.11 and also 

m 

i=0 
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by the induction hypothesis. Thus, 

and hence Taylor's Theorem with remainder yields 

+ !(■),£) - h^{-),e) = (B.IO) 

since and its derivatives are 0{1). This is the desired estimate of the first term 
in the right member of Eq. (]B.9|) . 

It remains to estimate the second term, G'(-, £:) in the right member of 
Eq. flRQj) . We recall that G= (^^j), where / and g are 0{1) in general. Hence, the 

first component of (?(■, hm+i{-),e) is plainly 0{e). The second component is as well, 
since Lemma [2TT] implies that /im+i,o = ^[o] and hence that g{-,hm^i{-),e) = 0{e), 
also. Therefore, 

G{;h^^i{-),e) = 0{e). (B.ll) 



Combining the estimates ( IB. 101) and (IB. Ill) , we see that the right member of 



Eq. (IB. 91) is which is the desired result. 

Finally, the solution of the condition Lm+i = yields an Ng— dimensional man- 
ifold Cm+i, as may be shown using the Implicit Function Theorem and [HI Theo- 
rem 1.13]. It suffices to show that 

det (DyLm+i) (-, hm+ii-)) ^ 0. 
Lemma [B.ll yields a leading order formula for DyL^+i, 

{DyL^^,) {z) = {-e-'H (Dyg)^)"'^' + O (e, ||^7o(^)||) ■ 

Here, z is a general point and {■)o{z) = {■){z,0). Next, we showed above that 
h{m+ifi) = ho. Recalling, then, Eq. (12.51) . we obtain 

(DyLrn+i) (x, hrn+lix)) = [-£-^H {Dyg)^]""^^ + for all x e K, 

where {Dyg)o = {Dyg){x, ho{x),0). Thus, 

det (DyLm+i) {x, hm+i{x)) ^ 0, for all x G K, 
by normal hyperbolicity and the proof is complete. I 
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